function [IRF_Gibbs, AA, non_stab] = IRF_BVARGibbs(gamma_g,Su_g,p,N,Gibbs,horizon);
%%
non_stab=0;
IRF_Gibbs=zeros(N,N,horizon,Gibbs);
AA=[];
J = zeros(N,N*p); J(:,1:N) = eye(N);
J=sparse(J);
for jg = 1:Gibbs
    A = zeros(N*p);
    A(1:N,:) = gamma_g(1:end-1,:,jg)';
    A(N+1:end,1:N*(p-1)) = eye(N*(p-1));
    vecModulus = abs(eig(A));
    AA(:,:,jg)=A;

    if max(vecModulus)>1
        non_stab=non_stab+1;
    end; %else
    B0 = chol(Su_g(:,:,jg))'; %% Lower triangular
    A=sparse(A);
    AJ=eye(size(A,1));
    for j = 0:horizon
        IRF_Gibbs(:,:,j+1,jg) = (J*AJ*J')*B0;
        AJ=AJ*A;
    end
    %end
end
